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Based directly on the microscopic lattice dynamics, a simple high temperature expansion can 
be devised for non-equilibrium steady states. We apply this technique to investigate the disordered 
phase and the phase diagram for a driven bilayer lattice gas at half filling. Our approximation 
captures the phases first observed in simulations, provides estimates for the transition lines, and 
allows us to compute signature observables of non-equilibrium dynamics, namely, particle and energy 
currents. Its focus on non-universal quantities offers a useful analytic complement to field-theoretic 
approaches. 
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I. INTRODUCTION 

Many-particle systems in a state of thermal equilibrium are the exception, rather than the rule. Physical reality 
is overwhelmingly in a far-from-equilibrium state. Examples range from living cells and weather patterns to ripples 
on water and sand. As we leave the framework of standard Gibbs ensemble theory for equilibrium systems, we have 
to search for new avenues and tools, seeking to understand and classify non-equilibrium behavior. As a first step 
along this road, the study of the simplest generalizations of equilibrium systems, i.e., non- equilibrium steady states 
(NESS), has been particularly fruitful [1,2]. Progress has relied predominantly on simulations, mean-field theory 
and renormalization group analyses for simple model systems. A class of models which exhibit especially interesting 
behavior are driven diffusive systems. Microscopically, these are lattice gases, consisting of one or more species 
of particles and holes, whose densities are conserved. An external driving force, combined with suitable boundary 
conditions, maintains a NESS. In the simplest case [3], a uniform bias, or drive, E, is imposed on an Ising lattice gas 
such that a nonzero steady-state mass current is induced. This model differs significantly from the usual Ising model: 
it displays generic long-range correlations [3-5], and belongs to a non-equilibrium universality class [6] with upper 
critical dimension dc — The ordered phase is phase-separated into two strips of high vs low density aligned with 
the bias. In contrast to equilibrium, bulk and interfacial properties are inextricably intertwined here [7]. 

To avoid the complications due to the presence of interfaces, a bilayer structure was suggested [8]: in the two- 
dimensional case, a second lattice was introduced, allowing for particle-hole exchanges between each site and its 
mirror image. This bilayer system is half filled with particles, and both layers are driven in the same direction. 
In the absence of any energetic couplings between the two layers, it was hoped that typical ordered configurations 
would show homogeneous densities on each layer, one almost full and the other nearly empty. Remarkably, however, 
this expectation proved too naive: Monte Carlo simulations [9] showed a sequence of two phase transitions, as the 
temperature is lowered: the first transition takes the system from a disordered (D) phase to a strip-like (S) structure 
showing phase-separation within each layer, with interfaces parallel to the drive and 'on top of one another. The 
anticipated "full-empty" (FE) phase, with uniform densities on both layers, only emerges after a second transition 
which occurs at a lower temperature. Once an interaction J, of either sign, between nearest neighbors on different 
layers is introduced, the full phase diagram in (J, E) space can be mapped out [10,11], using Monte Carlo simulations. 
As one might expect, the S (FE) phase dominates for attractive (repulsive) cross-layer coupling J. Remarkably, 
however, there is a small but finite region where the S-phase prevails even though the cross-layer coupling is weakly 
repulsive (cf. Fig. 1). The presence of this domain puts the two transitions, observed for J = 0, into perspective. We 
note for completeness that universal properties along the lines of continuous transitions have been analyzed in [12] 
with the help of renormalized field theory. 

To provide additional motivation for the study of layered structures, we note that multilayer models have a long 
history in equilibrium statistical mechanics [13-15]. On the theoretical side, they allow for the study of dimensional 
crossover [16]; on the more applied side, they provide natural models for the analysis of intercalated systems [17], 
interacting solid surfaces or thin films [18]. Since intercalated systems are often driven by chemical gradients or 
electric fields, to speed the diffusion of foreign atoms into the host material, it is quite natural to study driven layered 
structures. 
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Simulations relie, of course, on discrete lattice models. In contrast, field theories operate in the continuum, and 
thus, all discrete degrees of freedom have to be coarse-grained before these powerful techniques can be applied. In 
the process, non-universal information is lost, such as, e.g., the location of transition lines in the phase diagram. It is 
therefore desirable to identify a second analytic approach which is based directly on the microscopic model and thus 
complements both, simulations and continuum theories. Fortunately, high temperature expansion techniques [19] can 
be generalized to interacting driven lattice gases [4,20,21]. For the single-layer case, two-point correlation functions 
can be computed approximately [4,20] and display the expected power law decays in the steady state. With some 
care, the approximate location of order-disorder transitions can be extracted and compared to simulation results [20]. 
Given the nature of the approximation, quantitative accuracy cannot be expected, but the qualitative agreement of 
data and approximation is remarkably good. 

While the high temperature expansion is quite successful for the usual driven lattice gas, it is not clear to what 
extent it is capable of capturing the main features of other driven systems. This motivates the work presented in this 
paper, namely, the analysis of the bilayer system with this technique. Within a first-order approximation, we compute 
the two-point correlation functions and several related quantities, such as the particle current and the energy flux 
through the system. We extract the approximate location of the continuous transition lines and compare our results 
to the Monte Carlo data. As in the single-layer case, the qualitative features of the transition lines are reproduced as 
well as can be expected. Some limitations of the method will be discussed. 

This paper is organized as follows. We first introduce the bilayer model. After a brief summary of the high 
temperature expansion, we derive the closed set of equations satisfied by the two-point functions. We then obtain the 
solutions and extract the transition lines. Next, we show how the mass and energy currents through the system can 
be expressed in terms of pair correlations. We conclude with some comments and open questions. 



II. THE BILAYER MODEL 



X 2. Each lattice site f = {x,y,z), with x,y = 1,2,...,L and z = 0,1, 
wc also use lattice gas language, mapping spins into particles or holes, via 



A variant of the driven Ising model [3], the model consists of two square lattices, one stacked above the other, 
resulting in a bilayer structure of size 
carries a spin variable ,s(r) = ±1. Often 
s{r) = 2n{r) — 1. The local occupation variable n(f) takes the values 1 or 0, indicating whether a particle is present 
or not. The total magnetization, X]^s(r), is fixed at zero so that the Ising critical point can be accessed. Within each 
layer, nearest-neighbor spins interact through a ferromagnetic exchange coupling Jq > 0; in contrast, the cross-layer 
interaction J, which couples spins s{x,y,0) and s{x,y, 1), can take both signs. These choices are motivated by the 
physics of intercalated systems [17]. Thus, the Hamiltonian of the system can be written in the form 



H = - Jo^^ s{x, y, z)s{x', y' , 2) - 2 J ^ s{x, y, 0)s{x, y, 1) 



(1) 



where X^nn denotes the sum over all nearest-neighbor pairs (x, y, 0) and (x', y' , z) within the same plane. A heat bath 
at temperature T is coupled to the system, in order to model thermal fluctuations. We use fully periodic boundary 
conditions in all directions; hence the factor of 2 in front of the cross-layer coupling J. 

In the absence of the drive, particles hop to empty nearest-neighbor sites according to the usual Metropolis [22] rates, 
min {1, cxp (— /3A_ff)}, where AiJ is the energy difference due to the jump. Respecting the conservation of density, 
the phase diagram of this system is easily found. At high temperatures, a disordered phase persists, characterized by 
correlations which fall off exponentially. At a critical temperature Tc(J), a continuous transition occurs into the S 
(FE) phase for J > (J < 0). At J = 0, the critical temperature takes the Onsager value [23] rc(0) = 2.269... Jq/Zcb. 
For finite J, Tc{J) is even in J, due to a simple gauge symmetry, and increases monotonically with | J|. For J —>■ 
±00, nearest-neighbor spin pairs, with the partners located on different layers, combine into dimca's who couple to 
neighboring dimers with strength 2Jo. As a result, the critical temperature approaches the limit Tc(±oo) = 2Tc(0). 
The line J = 0, T < Tc (0) is a line of first-order transitions between the S and FE phases. It ends in a bicritical point 
at J = 0, T = Tc{0). 

To drive the system out of equilibrium, we apply a bias (an "electric" field) E along the positive x-axis. The 
contents of two sites, r and f+ a, separated by a (unit) lattice vector o, are exchanged according to the rate 



c(f, r + a; {s}) = min |l,exp —PAH + f3a- E (n(f) — n{r + a))j | 



(2) 



The argument {s} reminds us that the rate depends on a local neighborhood of the central pair. Due to E, particle hops 
against the drive become unfavorable. In conjunction with periodic boundary conditions in the x- and y-directions, 
the system settles into a non-equilibrium steady state with a net particle current. 
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The phase diagram, resuhing from Monte Carlo simulations at Jo = 1 and infinite E, is shown in Fig. 1. The 
same phases and transitions are found, but the bicritical point and its attached first order line are shifted to higher 
T and into the J < region. Thus, the S phase is observed to be stable in a finite window of negative interlayer 
coupling, so that two transitions must occur along the J = axis. This discovery represents the most unexpected 
new characteristic of this driven diffusive system. We also note the decrease of the critical temperatures for very large 
\J\. In a recent paper [11], this phase diagram was extended to include unequal intra-layer attractive couplings. In 
this case, the bicritical point is shifted even further into the negative region of J as the coupling transverse to the 
bias increases. 

We now turn to the analysis of this model in terms of a high temperature expansion. 
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FIG. 1. The phase diagram of the driven bilayer lattice gas from Monte Carlo simulations [10], at Jo = 1 and infinite E. T is 
measured in units of the Onsager temperature, T'c(O). The domains of the three phases are indicated. The inset demonstrates 
the shift of the bicritical point to negative J. Open (filled) circles indicate continuous (first order) transitions. 



III. HIGH TEMPERATURE EXPANSION 



The dynamics underlying the Monte Carlo simulations is easily expressed via a master equation. The latter provides 
a convenient starting point for a high temperature expansion. For simplicity, we take the thermodynamic limit within 
each plane, i.e., L — > oo. Following [4], we first derive the equations of motion for the two-point functions. By virtue 
of the familiar hierarchy, they are coupled to the three-point functions; however, we will argue that these are negligible 
(while non-zero, they are numerically rather small), so as to arrive at a closed system of equations for the two-point 
correlations. Temperature appears in these equations through the rates, via the combinations 0J, /3Jo, and (3E. To 
preserve the non-equilibrium nature of our dynamics, we expand in /3J and /3Jo, keeping PE finite. Technically, this 
requires that E always dominates the energetic contribution, i.e., E > AH for all jumps along E. To first order, a 
linear, inhomogeneous system of equations results, which can be solved exactly [20] and forms the basis of our analysis. 



A. The equations of motion and their solution. 



Before turning to any detailed calculations, let us introduce the key quantities. The two-point correlation function 
is defined as: 

G{f-f') = {s{r^s{r')) (3) 

where (•) denotes the configurational average. Due to translation invariance, G depends only on the difference of the 
two vectors. Moreover, G is invariant under reflection of one or several lattice directions; e.g., G{x, y, z) = G{—x, y, z), 
etc. The correlation function at the origin is obviously unity, G(0) = {s^{r)) = 1- We also introduce the Fourier 
transform of G, i.e., the structure factor. 



oo 



(4) 



Since wc take the thermodynamic Umit L ^ oo, the wave vectors k and p are continuous, but restricted to the first 
BriUouin zone [— tt, tt], while q is discrete, taking only the two values and tt. For completeness, we also give the 
inverse transform, 



= j 5'(A;,p,Q)e*('=^+f^+«^) (5) 

where the second line just defines some simplified notation. 

To set up the high temperature expansion, we first define the actual expansion parameters of our theory, namely 

= /3 Jo 

K = f3J (6) 

For K — Kf) = 0, the steady-state distribution is exactly known [24] to be uniform for all E: P* oc 1, so that we 
are expanding about a well-defined zeroth order solution. The correlation functions and structure factors for this 
limit are trivial, namely, G{r) = where 5 denotes the Kronecker symbol, and S{k,p,q) = 1. Returning to the 

interacting case, we note that G{f), for r 7^ 0, is already of first order in the small parameter. Similarly, we can write 
the structure factor as a sum of two terms. The first term is just the zeroth order solution, while the second, S, carries 
the information about the interactions, 

S{k,p,q) = l + Sik,p,q) (7) 
so that we can recast G(f), for rj^O, in the form 

G{x,y,z) = j S{k, p, q)e'ik-+Pv+i-) ioxx,y,zi^Q (8) 

The exact equations of motion for G are easily derived from the master equation [4] : 

I {s{f)s{r ')) = ^ {s{f)s{r ') [six)six ')-l]c {x, x {s})) (9) 

x,x' 

Here, the sum runs over nearest-neighbor pairs (.x, x ') such that x E {r, r '} but x ' ^ {r, r '}. Stationary correlations 
are obtained by setting the left hand side to zero. Clearly, jumps along and against all three lattice directions will 
contribute to the right hand side of Eq. (9). 

To proceed, let us write the jump rates in a form which makes their dependence on the spin configuration {s} 
explicit, so that the configurational averages in Eq. (9) can be performed. For infinite drive, a particle jumps along 
the field with rate unity, but never against it, so that the transition rates parallel to the field can be written as: 

cf {f,r + x;{s}) = ^[s{r)-s{f+x) + 2] (10) 

Here, i is a unit vector in the positive x-direction. In the case of finite drive, our restriction E > AH ensures that 
jumps along E still occur with unit rate, while those against E are suppressed by a factor of exp [— /? (Aif + 
Defining 

e = e-^^ (11) 

Eq. (10) must be amended to 

C|| {r,f+ x; {s}) = ^ [s{r) - s{f+ x) + 2] + ^ [s{f+ x) - s{f} + 2] exp{-(3AH) (12) 

Transverse to the field we have two jump rates, corresponding to the two transverse directions {y and z). Both of 
these are regulated by the energy difference due to a jump: 

c±{f,f+a;{s}) = min{l,exp(-/3Aif)} (13) 
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We are now ready to expand the rates in powers of K and Kq while keeping e finite: 



C|, if, f+ x; {s}) = i [s{f} - s{r + x) + 2] + ^ [s{f+ x) - s{f) + 2] (1 - /3AiJ) + 0{p'') (14) 

(f, r + a; {s}) = 1 + /3c2 (f, r + a; {s}) + 0{l3^) (15) 

with 

C2{r,f+a;{s}) = ~{AH+\AH\) (16) 

Given these simple forms for the rates, we can now derive the equations of motion satisfied by the pair correlations 

directly from Eq. (9), following [4]. A few details arc outlined in the Appendix. Keeping only corrections to first 
order in K, Kq and neglecting three-point correlations, wc obtain a closed set of linear equations for G{x, y, z). Since 
the dynamics is restricted to nearest-neighbor processes, it is not surprising that the equations involve an anisotropic 
lattice Laplacian acting on G{x,y,z). For x,y,z near the origin, the Laplacian may include the origin and will thus 
generate inhomogeneities in the system of equations. The detailed form depends on the chosen boundary conditions, 
and, of course, on the three parameters K, Kq, and e. Below, wc; show the set of equations for fully periodic boundary 
conditions. The first three equations result from nearest neighbors of the origin, r = (1. 0, 0), (0, 1, 0), and (0, 0, 1): 

dtG{l, 0, 0) = (1 e)[G{2, 0, 0) - G{1, 0, 0)] + 4[G'(1, 1, 0) - G(l, 0, 0)] 

+4[G(1, 0, 1) - G(l, 0, 0)] + 2eKo + 8Ko 
dtG{0, 1, 0) = 2(1 + £)[G(1, 1, 0) - G(0, 1, 0)] + 2[G(0, 2, 0) - G(0, 1, 0)] 

+4[G(0, 1, 1) - G(0, 1, 0)] + 4sKq + 6Ko (17) 
dtG{0, 0, 1) = 2(1 + £)[G(1, 0, 1) - G(0, 0, 1)] + 4[G(0, 1, 1) - G(0, 0, 1)] 

+8K + 8eK 

By virtue of invariance under reflections, these equations also hold for the other nearest neighbors r = (—1,0,0), 
(0,-1,0), and (0,0,-1). The following three equations arise from the next-nearest neighbor sites, r = (1,1,0), 
(0, 1, 1), and (1, 0, 1), and their reflections: 

dtG{l, 1, 0) = (1 + e)[G{2, 1, 0) + G(0, 1, 0) - 2G(1, 1, 0)] + 2[G(1, 2, 0) + G(l, 0, 0) 

-2G(1, 1, 0)] + 4[G(1, 1, 1) - G(l, 1, 0)] - 2Ko - 2eKo 
dtG{0, 1, 1) = 2(1 + £)[G(1, 1, 1) - G(0, 1, 1)] + 2[G(0, 2, 1) + G(0, 0, 1) - 2G(0, 1, 1)] 

-h4[G(0, 1, 0) - G(0, 1, 1)] - 4[Kq + K] (18) 
5tG(l, 0, 1) = (1 + £)[G(2, 0, 1) + G(0, 0, 1) - 2G(1, 0, 1)] + 4[G(1, 1, 1) - G(l, 0, 1)] 
-h4[G(l, 0, 0) - G(l, 0, 1)] - AeK - 4Ko 

Increasing the separation of the participating sites further, to f* = (2, 0, 0) and (0, 2, 0), wc obtain: 

dtG{2, 0, 0) = (1 + e)[G(3, 0, 0) + G(l, 0, 0) - 2G(2, 0, 0)] + 4[G(2, 1, 0) - G(2, 0, 0)] 

+4[G(2,0,l)-G(2,0,0)]-2eifo 
dtG{0, 2, 0) = 2(1 + e)[G(l, 2, 0) - G(0, 2, 0)] + 2[G(0, 3, 0) + G(0, 1, 0) - 2G(0, 2, 0)] (19) 

+4[G(0,2,1)-G(0,2,0)] -2ii:o 

And finally, all G's with |a;| + \y\ + \z\ > 2 satisfy homogeneous equations: 

dtGii,j, k) = il + e) [G{i +l,j,k) + G{i -l,j,k)- 2G(i, j, k)] 

+2[G{i,j + l,k) + G{i,j-l,k)- 2G{i,j, k)] + 4[G{i, j, k - 1) - G{i,j, k)] (20) 

The last equation contains the full anisotropic lattice Laplacian, acting on G{i,j,k), without any inhomogeneities 
being generated. Wc note, for further reference, that the right hand sides of Eqns (17-20) contain contributions from 
exchanges along and against the three lattice directions. Starting from Eq. (9), it is of course easy to keep track of 
terms originating in transverse vs parallel jumps. Below, this distinction will become important when we turn to 
energy currents. 

To solve this system, we closely follow the method presented in [20]. Returning to Eq. (7), we need to focus only 
on S, since this quantity carries the information about the interactions. Recalling Eq. (8), we first express G through 
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its Fourier transform S, exploiting translation invariance and linearity. Then, we invoke the completeness of complex 
exponentials to project out an equivalent set of (algebraic) equations for S. 

To follow through with this program, we first define the anisotropic lattice Laplacian in Fourier space, 

5{k,p,q) = 2(1 + £)(!- cos 4(1 -cosp)+ 4(1- cos g) (21) 

and second, introduce the three (as yet unknown) quantities: 

ji = y ^(1 -cosfc) 

l2= j 5(1 - cosp) (22) 
l3= J 5(1- cos g) 

With these definitions, the system can be expressed in terms of S, resulting in: 

2eKo + 8Ko = J S5ex-p{ik) + (1 + e)/i 
AeKo + 6Ko = j S6exp{ip) + 2h 
%eK + %K = j S5 exp{iq) + 4/3 
-2eKo-2Ko = j S6exp{i{k + p)) 
-4:{Ko + K) = J S5 exp(i(p + q)) (23) 
-4eK-4Ko = J S6exp{i{k + q)) 
-2eKo = j SSexp{2ik) 
-2Ko = j S6exp{2ip) 

= y S6 ex.p{i{kx + py + qz)) foTc |x| + |y| + |^| > 2 

To proceed, we treat /i , /2 , ^3 for the time being as simple coefficients and move them to the left-hand side. Finally, 
we need one additional equation iov x = y = z = 0, which is easily obtained: 

j SS = j 5[2(l + £)(l-cosfc)+4(l-cosp)+4(l-cosg)] = 2(1 + e)7i + 4/3 + 4/3 

Now, we are ready to invoke the completeness relation for complex exponentials, namely, 

exp[i{kx +py + qz)] = 2{2TTfS{k)S{p)Sg,o (24) 



x,y,z 

which allows us to solve for S: 



where 



Qti, ^ L{k,p,q) 

S{k,p,q) = — (25) 

S{k,p,q) 



L{k,p,q) = 2(1 + £) (1 - cosfc)/i + 4 (1 - cosp) /a + 4 (1 - cosg)/3 

+ i2eKo + 8K0) 2 cos fc + (4£i^o + Si^o) 2 cosp 

+ {8sK + 8K) cos q - {2eKo + 2Kq) 4 cos k cosp (26) 
- {AeK + AKq) 2 cos A; cos g - A{Ko + K)2 cosp cos q 
-AeKq cos 2k - AKq cos 2p 
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However, Eq. (25) is not yet a fully explicit solution for S, due to the appearance of the three integrals Ii, I2, and I3 
in L. To determine these three coefficients, we need three linearly independent equations. One of these equations is 
given by the value of G at the origin, 1 = G{0, 0, 0) = J{1 + S), and the remaining two can be obtained directly from 
the definitions of Ji and I3 in Eq. (22): 







L{k,p, q) 



d{k,p,q) 

= -/s + / 777 r 1 - cosg) 

J S{k,p,q) 

After inserting Eq. (26) for L, this leads to a set of three inhomogeneous, linear equations for the three unknowns h, 
I2, and I3, which are easily solved. Since the details of the associated matrix inversion are straightforward but tedious, 
we relegate a few details to the Appendix. We just note the following overall features: (z) All three coefficients are 
functions of K, Kq, and e; {ii) for the whole range of fields e and for Kq = 1 and K = ±1 (attractive and repulsive 
inter-layer interactions), h and I2 are negative, while I3 is positive for ii' = — 1 and negative for K = +1. 
This concludes the calculation of the structure factor. To summarize, we find 

S{k,p, q) = l + + 0{K^Kl KKo) (28) 

d{k,p, q) 

Even at the lowest nontrivial order, this solution carries a significant amount of information about the phase diagram 
of our system. In particular, we can extract an approximate shape of the critical lines, as we will show in the following. 



B. The critical lines. 



The location of a continuous phase transition is marked by the divergence of a suitably chosen structure factor, 
as a function of the external control parameters. For example, we can locate the order-disorder phase transition of 
the usual, two-dimensional Ising model by seeking those values of temperature (at zero magnetic field) for which the 
structure factor, S{k), diverges. In the absence of a conservation law, the only singularity occurs at the Onsager 
temperature if = 0, indicating that the system orders into a spatially homogeneous state. For a lattice gas, however, 
S{Q) is fixed by the conservation law, and we need to seek the onset of phase separation, i.e., the emergence of 
macroscopic spatial inhomogeneities in the system. In this case, singular behavior occurs in limg^QS'(fc), provided 
the system is half-filled and tuned to the Onsager temperature. 

For the bilayer system, we need to locate, and distinguish, two types of continuous transitions, namely, from 
disorder (D) into the strip (S) and the full-empty (FE) phases, respectively. Since the D-S transition is marked by 
the appearance of phase-separated strips in each layer, aligned with the driving force, it can be located by seeking 
singularities in limp^o S{0,p, 0). In contrast, the D-FE transition exhibits homogeneous, but opposite magnetizations 
in the two planes, so that it can be found by considering S{0, 0, tt). In fact, these two structure factors were precisely 
the order parameters chosen in the MC studies [10]. 

Yet, another subtlety must be considered: in a typical high temperature expansion such as ours, only a finite 
number of terms can be computed. Hence, any perturbative result for the structure factor must be finite, and instead, 
the radius of convergence of the expansion must be estimated. Even this is not practical here, since we have only two 
terms of the series. To circumvent these restrictions [20], we extract the singularity by looking for zeros oi S~^, to 
first order in K and Kq. 

Starting from Eq. (28), we obtain 

S-\k,p, q) = l- + 0{K\Kl KKo) (29) 

d{k,p,q) 

and seek to locate the zeros of limp^o S~^{0,p, 0) for the D-S transition, and of 5~"'^(0, 0, tt) for the D-FE transition. 
Of course, wc should ensure that these are the first zeros which are encountered upon lowering the temperature. 
Therefore, we consider, more generally, the behavior of S~^{k,p,q) at small k,p and fixed q. The denominator of 
Eq. (29), being the lattice Laplacian, is positive definite: 

lim 6{k,p,q)= (1 -|- e)/;;^ -|- 2/ + 4(1 - cosg) 0(A;*,/, fcV) (30) 
fe,p— >o 
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and vanishes only at the origin. Similarly, we obtain 

lim L{k,p,q) = 16Ko{l- cos q)+ A {1- cos q)l3 (31) 
fe,p— »o 

[(1 + s)Ii + IOsKq - 4Kq + [AeK + AKq) cos q] 
V [2/2 + + {AKq + AK) cosg] + 0(fc^/, fcV) 

We note, briefly, that the anisotropic momentum dependence of numerator and denominator leads to power law 
correlations in the x- and y-directions [4,20,25]. Combining our results so far, it is apparent that the zeros of are 
identical to those of (5 — L in Eq. (29). To simplify notation, wc write 

lim [5{k,p,q) - L{k,p,q)] = ^{q)k'' + 2T^{q)p^ + Atz{1 - cosg) + 0(fc^p^ fcV) (32) 
fe,p— »o 

and read off 

T|| (g) = (1 + e) (1 - 7i) - lOeKo + AKq - {AeK + AKq) cos q 

T± {q) = l-h- 3Ko -{2Ko + 2K) cos q (33) 
= 1 - 73 - 47^0 

In a field-theoretic context [12], these quantities play the role of diffusion coefficients: T|| and t± control the in-plane 
diffusion in the parallel and transverse directions, respectively, while controls the cross-plane hopping. 

For high temperatures, i.e., small values of Kq = (3 Jo and K = (3 J, all three r-coefficients are positive. Seeking 

zeros of these expressions, as Kq and K increase, wc need to consider the two cases q = and q = ir separately. For 
g = 0, we find that tj_(0) has a single zero at a critical /?f , for given Jq, J and e. At these parameter values, r|| (0) and 
Tz remain positive. Similarly, for q = w, the coefficient Tz is the one which vanishes first as /3 increases, reaching zero 
at a critical /3f ^. Converting into temperatures, we obtain two functions, T^f ( Jo, J, e) and Tf^{Jo, J, e), and we need 
to identify the larger of the two: If max [T^^( Jq, J, e),T^^{Jo, J, £)] = T^{Jo, J, e), the order-disorder transition is of 
the D-S type. Otherwise, if max [T^( Jq, J, e), T^^( Jo, J, e)] = Tf^{Jo, J, s), the FE phase is selected upon crossing 
criticality. 

While the two critical lines, T^{Jo, J, e) and T^^{Jq, J,s), can in principle be found analytically, the details are not 
particularly illuminating. Instead, we present a range of numerical results below. For example, for infinite E (e = 0), 
we obtain 

kBT^{Jo, J, 0) = 4.39./0 + 2.11 J (34) 
kBT^^{Jo, J, 0) = 4.14Jo - 1.36J 

For finite E with e = exp(— /37J) = 0.5, all coefficients decrease and we find 

ksT^i Jo, ^,0.5) = 4.I5J0 + 2.03J 
fcsTf ^(Jo, J, 0.5) = 4.O5J0 - 1.70J 

In each case, the bicritical point is defined through the solution of T/(Jo, J, e) 
also quote the equilibrium {E = 0) results (see the Appendix for details): 

fcBTf(Jo, J,l) =4Jo + 2J 
A;Brf^(Jo, J,l) =4Jo-2J 

which exhibit the expected J — J symmetry. 

Recalling that the MC simulations were performed at fixed, positive in-plane coupling Jo, wc need to consider only 
the dependence on the cross-plane coupling J which may take either sign. All of our results show that, for positive 
J, the D-S transition dominates while, for sufficiently negative J, a D-FE transition is found. In the following, we 
discuss the non-equilibrium {E ^ 0) results in more detail. 

Fig. 2 shows the critical lines for two typical values of the parameters, e = 0.5 and Jo = 1. Being the result 
of a first order approximation, the critical lines must of course be linear in J. Therefore, quantitative agreement 
with the simulation data cannot be expected; nevertheless, several important features are reproduced: the existence 
of two ordered phases and the shift of the bicritical point to higher values of T and negative J. As a result, the S 
phase survives for small, negative J, despite being energetically unfavorable. This phenomenon can be explained 
qualitatively [10] by noting that long-range negative correlations transverse to E dominate the ordering process for 



(35) 

= T^^{Jq, J,e). For comparison, we 

(36) 
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positive J, and this mechanism continues to be effective for a small region of negative J. For large and negative J, the 
disordered state orders into the full-empty (FE) phase, characterized by the planes being mainly full or empty. Finally, 
we comment on the dependence of the critical lines, specifically T^{1, l,e) and Tf^{l, —l,e), on the field parameter 
e = exp {—pS), shown in Fig. 3. For E = 0, both temperatures are equal, by virtue of the J — » — J symmetry of 
the equilibrium system. As the field becomes stronger, the critical temperature of the D-S transition increases, in 
contrast to the critical temperature of the D-FE transition which decreases. This behavior agrees qualitatively with 
the trend observed in the simulations [10,26]. 




-0.1 



-0.05 



0.05 



0.1 



FIG. 2. The dependence of the critical temperature, max [rf ( Jo, J, e), ^( Jo, J, e)] , on the cross-layer coupling J, for 
e = 0.5 and Jo = 1. Filled diamonds (open squares) indicate that the transition is of the D-S (D-FE) type. Critical temperatures 
for J > 0.1 (J < —0.1) can be obtained by linear extrapolation of the D-S (D-FE) branch. 

There are several other quantities of physical interest which are immediately related to the two-point correlations, 
such as the steady-state particle and energy currents. To extract these, we first discuss the inverse Fourier transform 
of the structure factor, focusing specifically on the nearest-neighbor correlations. 
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FIG. 3. The critical temperatures, T^(l, l,e) corresponding to a D-S transition (filled diamonds) and Tf^(l,- 
sponding to a D-FE transition (open squares) as functions of the field parameter e = exp{—f3E). 



-l,e) corre- 



C. Related physical observables. 

Nearest-neighbor correlations. These are easily found from our solution for the structure factor, Eq. (28). For 
example, the nearest-neighbor correlation in the field direction is given by: 



9 



G(1,0,0) = J S{k,p,q)cosk + 0{K^,K^,KKo) (37) 
Since J S = hy virtue of G(0, 0, 0) = 1, we obtain 

G{1,0,0) = -J S{k,p,q){l-cosk)+0{K^,KlKKo) = -Ii (38) 

and similarly, 

G(0,l,0) = -72 (39) 
G(0,0,l) = -73 

These three integrals are already known since they were required for the discussion of the critical lines. Specifically, 
for e = 0.5 we find, neglecting corrections of 0{K^, Kq.KKq): 

G(l, 0, 0) = 0.949i^o + O.OBOis: 

G(0, 1, 0) = 0.849i^o - QMAK (40) 
G(0,0, 1) = -0.055ii'o + 1.702ii' 

For reference, we also quote the first order results for the equilibrium {E = 0) correlations: 

G'^«(l,0,0) = G^«(0,l,0) = ii'o (41) 

In the following graphs (Fig. 4a-c), we show the drive dependence of the three nearest- nez(7/i6or correlation functions, 
at i^o = 1 and K = ±1, to illustrate their behavior in two typical domains (attractive and repulsive cross-layer 
coupling). Of course, these values of K and Kq are not "small" , but in a linear approximation they just serve to fix a 
scale. Consistent with the interpretation of the drive as an additional noise which tends to break bonds, all correlations 
are reduced compared to their equilibrium value. Further, as the field is switched on, the J — » — J symmetry of the 
equilibrium system is broken, and the correlations for repulsive and attractive cross-layer coupling differ from one 
another. The details of how they differ provides some insight into the ordered phases which will eventually emerge. 

The first plot (Fig. 4a) shows the correlation function for a nearest-neighbor bond within a given plane, aligned with 
the drive direction. It is interesting to note that the correlations for repulsive cross-layer coupling are more strongly 
suppressed than their counterparts for attractive J. This feature becomes more transparent when we consider nearest- 
neighbor correlations transverse to the drive, but still within the same plane (Fig. 4b). For attractive cross-layer 
couphng, we note that G(l, 0, 0) is considerably enhanced over G(0, 1, 0), while the two correlations are roughly equal 
in the repulsive case. This indicates a tendency to form droplets of correlated spins which are elongated in the field 
direction for J = -1-1 while remaining approximately isotropic for J = — 1, hinting at the nature of the associated 
ordered phases (strip-like vs miiform within each layer). This picture is completed when we consider the cross-plane 
correlations G(0, 0, 1) (Fig. 4c): These are positive in the attractive, and negative in the repulsive case, demonstrating 
the tendency towards equal vs opposite local magnetizations on the two layers. Given that we have performed only 
a first order calculation, the results really carry a remarkable amount of information about the system. Encouraged 
by these observations, we now consider two other quantities, namely, the particle and energy currents. 
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FIG. 4. Nearest- neighbor pair correlations along the x- (a), y- (b), and z- (c) axes, as functions of the drive parameter e, 
for Kq = \K\ = 1. Filled diamonds (open squares) refer to ferromagnetic (antiferromagnetic) cross-layer coupling K = +1 
{K = -l). 
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The particle current. Due to the bias in conjunction with periodic boundary conditions, the bilayer system carries a 
net particle current, (j). Since only nearest-neighbor exchanges are possible, this current is proportional to the density 
(number per site) of available particle-hole pairs in the field direction. The transition rate C|| along this direction, 
given in Eq. (12), then counts the fraction of these pairs which will actually exchange per unit time. Specifically, in 
configuration {s}, the particle current can be written as 



1 ^ g(r) 



s{f+ x) 



2^2 



cJf,f+x;{s}) 



(42) 



For infinite E, this expression simplifies considerably, since jumps against the field will be completely suppressed. 

After a few straightforward algebraic manipulations, the average current can be expressed through the pair corre- 
lations along the field direction. To first order in K and Kq, we obtain 



1 



(43) 



which shows that it is non-zero only if £^ 7^ 0. Further, it takes its maximum value at infinite temperature and is 
reduced by (attractive) nearest-neighbor interactions. The graph (Fig. 5) shows the field-dependence of this current, 

for Kq = 1 and K = ±1. Since nearest- neighbor correlations along the field are much larger for positive J, indicating 
a predominance of particle-particle or hole-hole pairs, the current is reduced relative to the repulsive case. 
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FIG. 5. The average particle current, (j), vs the drive parameter e, for Kq = \K\ = 1. Filled diamonds (open squares) refer 
to ferromagnetic (antiferromagnetic) cross-layer coupling. 

Energy currents. Another interesting quantity associated with driven dynamics is the change in configurational 
energy during one Monte Carlo step. In the steady state, by definition, the average configurational energy is of course 
constant. However, particle- hole exchanges parallel to the field direction tend to increase the energy, since the drive 
can easily break bonds. In contrast, exchanges transverse to E are purely energetically driven and hence, prefer to 
satisfy bonds so that the energy decreases [1]. In summary, we have 



/dH\ / dH\ , , 

Even if a particle current were absent, the presence of energy currents would signal the non- equilibrium steady state. 

Since the configurational energy involves only nearest-neighbor bonds, it is obvious that only the time evolution of 
nearest-neighbor correlations plays a role in these two fluxes. Specifically, we have 

^"^ ( f ) „ = (I) „ °' 0) + 1' - (I) „ °' 1) 

where the subscript on the time derivatives reminds us to select only those processes which are due to parallel 
exchanges alone. These can be easily identified from the terms contributing to Eq. (9) or (17). Of course, there is 
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an analogous equation for (dH/dt) j_. Collecting the relevant contributions and multiplying both sides by the inverse 
temperature /3 to express everything in terms of Kq and K, we find: 

-2K {2(1 + £) [G(l, 0, 1) - G(0, 0, 1)] + 8sK} (46) 

The correlation functions spanning next- and next-next nearest neighbors which appear here can again be determined 

from our solution for the structure factors (see Appendix). The result, at Kq = 1 and if = ±1. is shown in Fig. 6 as a 
function of e. As expected, this flux is always non-negative and monotonically increasing as a function of E. We note 
that the current for X = — 1 is slightly larger than its counterpart for if = -1-1. Since it is a complicated function of 
the couplings and several correlations, we cannot offer a simple intuitive explanation of this property. 
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FIG. 6. The average energy current, L ^ {df3H/dt)^^, vs the drive parameter e, for Kg 
squares) refer to ferromagnetic (antiferromagnetic) cross-layer coupling. 
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IV. CONCLUDING REMARKS 



Based directly on the microscopic lattice dynamics, the high temperature series provides us with a simple analytic 
tool which complements field theoretic approaches. Even in a first order approximation, it is remarkable how many 
features of the MC results are at least qualitatively - reproduced. To summarize our results briefly, we derive, and 
solve, a set of equations for the stationary pair correlation functions and their Fourier transforms, the equal-time 
structure factors. By matching the series expansion of the latter with the expected critical singularity, we find two 
critical lines, separating the disordered phase from the strip phase (S) and the full-empty phase (FE), respectively. 
We also observe the shift of the bicritical point which marks the juncture of these two lines, in very good qualitative 
agreement with the simulations. To illustrate the non-equilibrium character of the steady state, we compute the 
particle current and the energy flux through the system. The particle current is determined by the nearest-neighbor 
correlations in the field direction, and takes its maximum value in the absence of interactions. Our findings for the 
energy current confirm intuitive expectations: parallel exchanges tend to lower, while transverse exchanges tend to 
increase, the configurational energy. 

A brief comment on boundary conditions is in order. Even though it is quite natural to use periodic boundary 
conditions in all lattice directions, it is not unreasonable to consider other choices, especially in the z-direction. To 
recall, periodicity in z implies that the site {x,y,0) is connected to the site {x,y, 1) via two bonds which enter into 
both the energetics and and the dynamics (i.e., there are two channels for a particle to move from one layer to the 
other). Alternately, we can choose open boundary conditions in z and consider only a single energetic bond and 
single dynamical channel between these two sites. Mixtures of these two cases can also be constructed: i.e., imposing 
periodic boundary conditions on the energetics, but allowing only a single channel for particle moves, or vice versa. 
The first (second) "mixed" case is reducible to the case of open (periodic) boundary conditions, with J replaced by 2 J 
(J/2). Even though details are not presented here, we did, in fact, compute the critical lines for different cross-plane 
boundary conditions. The main conclusions of our study, namely, the existence of the two continuous phase transitions 
and the shift of the bicritical point, hold for all of these variations. 
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The high temperature expansion presented here has two shortcomings. First, our resuhs provide no insight into 
the first-order transitions between the FE and S phases which were observed in the simulations. As in aU high 
temperature series, the first singularity which is encountered as T is lowered determines the radius of convergence. 
A low-temperature approach would be necessary to capture transitions between ordered phases. Second, our series is 
currently limited to just one nontrivial term. In order to compute the second order correction to the pair correlations, 
we would need to know the full stationary distribution, P*, to first order. Writing the stationary master equation in 
the form = LP* where L is the linear operator ( "Liouvillean" ) defined by the transition rates, this requires the full 
inverse of L, to zeroth order. Finding this inverse is a highly nontrivial (and as yet unsolved) problem. 

Inspitc of these drawbacks, the high temperature expansion is one of the few analytic tools which provide insight 
into non-equilibrium steady states. It is conceptually and mathematically straightforward, and - at least at the 
qualitative level - surprisingly reliable. Since it is based directly on the microscopic lattice dynamics, it still carries 
information about nonuniversal properties which would be lost upon taking a continuum limit. It is therefore a 
valuable complement to both simulations and field theoretic methods. 

Acknowledgements. We thank U.C. Tauber and R.K.P. Zia for fruitful discussions. Financial support from the 
NSF through the Division of Materials Research is gratefully acknowledged. 

APPENDIX A 

1. Equations for the two-point correlations. 

To illustrate the general procedure, we provide a few details here [4]. As an example, we choose the two-point 
correlation G(l, 1,0). We start with the equation of motion for the pair correlations: 

I {s{r)sr)) = E {simn [s{x)s{x') - 1] c {X, x'; {s})) 

where the sum nms over nearest-neighbor pairs (a;, x') such that x G {r, r '} but x' ^ {r, r '}. To obtain the equation 
satisfied by G(l,1.0), we choose, e.g. r= (0,0,0) and r' = (1,1,0). The two participating spins have a total of 8 
distinct nearest neighbors: 6 of these lie in the z = plane, and 2 are found on the z = 1 plane. One possible (x, x') 
pair is, for example, the pair x = r and x' = (1,0,0). The corresponding exchange occurs along the field direction, 
and hence, we must use the rate cy from Eq. (12). Considering only the contribution due to this particular pair of 
sites, we obtain 

dtGil, 1, 0) = i (.s(0, 0, 0)s(l, 1, 0) [s(0, 0, 0)s(l, 0, 0) - 1] 

X {[.s(0, 0, 0) - .s(l, 0, 0) + 2] + e [s(l, 0, 0) - s(0, 0, 0) + 2] exp(-/3Ai7)}) + ... 

Here, ... stands for the contributions due to all other possible {x,x') pairs which can be handled in an analogous 
manner. After multiplying out a few terms and neglecting 3-point functions, the expression above simplifies to 

atG(l,l,0) = i[G(0,l,0)-G(l,l,0)] 

+1 ([s(l, 0, 0)s(l, 1, 0) - s(0, 0, 0)s(l, 1, 0)] [s(l, 0, 0) - s(0, 0, 0) + 2] exp(-/3Aif )) + ... 

Note that we have replaced (s(l, 0, 0)s(l, 1, 0)) by the corresponding correlation function, G(0, 1, 0). Next, we expand 
exp(— /3AiJ) in powers of K and Kq, according to Eq. (14). The zeroth order contribution is easily accounted for, 
leaving us with the 0{P) correction: 

dtGil, 1, 0) = ^ (1 + e) [G(0, 1, 0) - G(l, 1, 0)] 

([s(l, 0, 0)s(l, 1, 0) - s(0, 0, 0)s(l, 1, 0)] [s(l, 0, 0) - s(0, 0, 0) + 2] (AH)) + ... 

The change in energy, AH, involves the nearest-neighbor spins of the selected pair. Since these terms are already 
explicitly first order in P, they are averaged using the zeroth order approximation to the steady state solution. The 
latter is uniform so that the averages are trivial. Collecting, we find that the contribution of this particular exchange 
to G(1,1,0) is 
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dtG{l, 1, 0) = i (1 + s) [G{0, 1, 0) - G{1, 1, 0)] - eKo + ... 

Including all the other (x,x') pairs, we arrive at the complete equation: 

dtG{l,l,0) = (1 + e)[G(2,l,0) + G(Oa,0) - 2G(1, 1, 0)] + 2[G(1, 2, 0) + G(l, 0, 0) 
-2G(1, 1, 0)] + 4[G(1, 1, 1) - G(l, 1, 0)] - 2Ka - IsKq 

Of course, this procedure is easily extended to any other two-point function. Moreover, it is straightforward to track 
which terms in the equations arise from parallel, and which from transverse, exchanges. This distinction is crucial for 
the computation of the energy fluxes. 

2. Matrix inversion. 

We seek to invert the system of equations (27) for the three unknown expressions Ji, l^-, and I^. We follow the 
method first outlined in [20]. With a little algebra, it becomes apparent that the coefficients and inhomogeneities in 
these equations involve integrals of the form 

/■ (1 — cos — cosp)™(l — cosg)" 



1= 1(1-' 

3 
2 



with non- negative integer /, m, n and I + m + n < 3. The task of determining these integrals is simplifled by a series 
of relations, namely, 

- = 2(1 + £)gioo + 4Qoio + 4Qooi 

(1 - cos k)- = 2(1 + £)Q2oo + 4giio + ^Qioi 

(1 - cosp)- = 2(1 + e)QnQ + 4Qo20 + 4Qoii 

(1 - cos = 2(1 + e)Q3oo + 4Q210 + 4Q201 

- = / (1 - cospf- = 2(1 + e)Qi2o + 4Qo30 + 4Qo2i 

(1 -cosfc)(l -cosp)- = 2(l + £)Q2io + 4Qi2o + 4Qiii 

(1 - cosg)- = 2(1 + e)Oioi + 4goii + 40oo2 

The computation of the remaining integrals, while tedious, is completely straightforward. Once these coefficients are 
known, Eqns (27) can be inverted. 

3. The equilibrium solution 

Since our expansion presumed E > AH, we may not simply set £ = 1 in our equations of motion for the two-point 
correlations. Instead, one should rederive the whole set carefully, noting the absence of the driving field. Of course, 
this is trivial, since only the first order term in the expansion of the equilibrium (Boltzmann) distribution is required 
here. In this order, only nearest-neighbor correlations can be nonzero, so that the only non-vanishing G's are 

G'^'^ (0,0,0) = 1 
G^«(±l, 0, 0) = G"«(0, ±1, 0) = Ko 
G^i{0,0,l) = 2K 

Performing the Fourier transform to structure factors and exploiting the boundary conditions in the 2;-direction, we 
find immediately that 
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2—0,1 :r,y— — oo 

= 1 + 2Kq (cos + cosp) + 2K cos g 

resulting in 

lim S-\k,p,q) = l-2Ko(2-\k^ -\p'\ ~ 2K cos q + 0{k\p^) 

/c,p— »0 \ 2 ^ ) 

= T{q)+0{k\p^) 

with 

T{q) = 1 - 4i('o - 2Kcosq 



For g = 0, this vanishes at 



and for q = it, the zero shifts to 



fcsTf = 4Jo + 2J 



fcsT<f ^ = 4 Jo - 2 J 



Thus, D-S transitions are observed for J > 0, and D-FE transitions dominate for J < 0. The bicritical point is located 
at ksTciJ = 0) = 4Jo. 

4. Other correlations near the origin. 

These are required to compute the energy fluxes along the parallel and transverse directions. Specifically, we need 
the following correlation functions: 1, 0), 0, 1), G(0, 1, 1), G{2, 0, 0), and G(0, 2, 0). We want to write these 
correlation functions in terms of the already calculated integrals Ii,l2, h and also in terms of the set of Qimn integrals 
defined earlier. 

We start with the definition for G{2, 0, 0) and substitute the expression for the structure factor: 

G(2,0,0)= J Scxp{2ik) = 2 J S'(l -cosfc)2 -4/i 

= 2 {2(1 + £)/iQ2oo + 4e(5Ko + 2if )Q30o 

+ 4(/2 + 5Ko + 2K)Q2w + 4 (/g + 4i^o) Q201 

- SKo{l + s)Q3io -S{K + Ko) Q211 

-8{eK + Ko) Q301 - 8eii'o(340o - &K0Q220 - 2/i} 

and similarly, 

G(0, 2, 0) = 2 J 8(1 - coapf - Ah 

= 2 {2(1 + e)/iQo2o + 4e(5i^o + 2K)Qi2o 

+ 4(/2 + 5Ko + 2K)Qo3o + 4 (/a + 4Ko) Q021 

- 8ii:o(l + e)Qi30 -SiK + Ko) Qosi 

-8{eK + Ko) Q121 - 8sKqQ22o - SKqQq^q - 2I2} 

The remaining correlation functions follow in the same way: 



G(1,1,0) = J S'(l-cosfc)(l-cosp)-/i -I2 



= {2(1 + e)/igiio + H^Ko + 2K)Q2io 

+4(/2 + 5Ko + 2X)gi2o + 4 (/g + 4Ko) Qm 

-8Ko{l + e)Q22o -8{K + Ko) Q121 

-8 {sK + Ko) Q211 - SeKoQsio - SJfoQiso - h - h} 
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G(1,0,1) = J S{1 - cosk){l - cosq) - h - I3 

= {2(1 + £)7iQioi + 4e{5KQ + 2K)Q2oi 

+4(12 + 57^0 + 2i^)Qiii + 4 (/g + AKq) Qw2 
-8ifo(l + e)Q2ii -8{K + Kq) Q112 
-8 [eK + Ko) Q202 - SeKoQaoi 
-8ifoQi2i - /i - /a} 

G(0, 1, 1) = 5(1 - cosp)(l - cos?) -h-h 

= {2(1 + e)/iQoii + 4£(5Jro + 2K)Qiii 

+4(/2 + S/To + 27r)Qo2i + 4 (/g + 4ii:o) Q012 
-8ifo(l + £)Qi2i - 8 (/C + i^o) Q022 
-8 {eK + i^o) Q112 - 8eifo02ii 
— Sii'oQosi — I2 — I3} 

After the additional Q-integrals have been determined, the energy currents are easily found. 
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